function f=ladder(Integrator)
% The Simulation of Ladder model of springs for Soft body dynamics
% Just call to execute.  
%Integrator values correspond to:
% 1 - Euler
% 2 - Runge-Kutta
% 3 - Verlet
% Settings for size shape, initial configuration and parameters can be modified from source

close all
clear all

% initialisations
N=10;%number of segments
n = 2*N;% number of point masses

m = 0.5; % mass of each point masses
kh =10*(2*2*pi)^2; % Stiffness of the horizontal springs
kv = 1*(2*2*pi)^2; % Stiffness of the vertical springs

lh =0.1*ones(1,N-1) ; %intrinsic length of horizontal springs
lv =  0.1./cumsum(ones(1,N)); %intrinsic length of vertical springs

tMax = 5; % max time of simulation in secs
deltaT = 0.5e-3; % step size of simulation in secs
b = 10; % kgm
initialAngle = 120/360*(2*pi);

%Setting just 10 runs of simulation (for debugging)
%tMax =1* deltaT;

%Function to set an initial curvature
[pV  pD ] = initialPosition(initialAngle , lh  , lv,N); 

%Storing the positions of first Ventral and Dorsal Masses
pDOrig =  pD(1,:); 
pVOrig = pV(1,:);

%p = [pV pD];

%velocities
pVDot = zeros(n/2,2);
pDDot = zeros(n/2,2);
pDot = [pVDot pDDot];
%pVDot(2,2)=-1;
%pDDot(1,1)=0.5;

%accelerations
pVDoubleDot = zeros(n/2,2);
pDDoubleDot = zeros(n/2,2);

%pDoubleDot = [pVDoubleDot pDDoubleDot];

%forces acting on each point mass
fD = zeros(n/2,2);
fV = zeros(n/2,2);

f = [fV fD];

% Dynamic Simulation by Euler Method

ctr = 0; 
% Fps
%fps=30;

i=2:N-1;	
while (ctr < tMax/deltaT)
%Measuring looptime with tic and toc
%	tic

% Plotting the points
    if (mod(ctr,1000*deltaT)==0) 
	
	%Plotting the position of the masses
	plot(pV(:,1),pV(:,2),'or',pD(:,1),pD(:,2),'ob')	

	%Drawing the connecting lines
	line([pV(1:(N-1),1)';   pV(2:N,1)'] ,[pV(1:(N-1),2)'; pV(2:N,2)'],'Color','r');
    	line([pD(1:(N-1),1)';   pD(2:N,1)'] ,[pD(1:(N-1),2)'; pD(2:N,2)']);
    	line([pD(1:N,1)';   pV(1:N,1)'] ,[pD(1:N,2)'; pV(1:N,2)'],'Color','g');
      
       % axis ( [-10 10 -10 10 ], "square");
         pause(0.01*deltaT);
         disp (ctr);
     end

	if(integator=1)
		%Integrating by Eulers Method
		%Forces Exterted by the springs
		[Ftd Ftv] = forceCalculator(pV,pD,pVDot,pDDot,N,kh,kv,lh,lv,b)	;
		
		%Calculating accelerations	
	        pDDoubleDot(1:N,:) = Ftd(1:N,:)./m;
		pVDoubleDot(1:N,:) = Ftv(1:N,:)./m;	
	
		%Calculating Positions
		pD(1:N,:) = pD(1:N,:) + deltaT*pDDot(1:N,:) +0.5* deltaT^2 *pDDoubleDot(1:N,:);
		pV(1:N,:) = pV(1:N,:) + deltaT*pVDot(1:N,:) +0.5* deltaT^2 *pVDoubleDot(1:N,:);
	
		%Calculating Velocitis
	       	pDDot = pDDot + deltaT *pDDoubleDot;
		pVDot = pVDot + deltaT *pVDoubleDot;		
    	end
        ctr = ctr+1;
        
        %% Fixing first Dorsal and Ventral Node to original Position
        pD(1,:)=pDOrig;
        pV(1,:)=pVOrig;

   %     toc
end